home *** CD-ROM | disk | FTP | other *** search
/ TeX 1995 July / TeX CD-ROM July 1995 (Disc 1)(Walnut Creek)(1995).ISO / graphics / gnuplot / standard.c < prev    next >
C/C++ Source or Header  |  1993-09-15  |  20KB  |  987 lines

  1. #ifndef lint
  2. static char *RCSid = "$Id: standard.c%v 3.50 1993/07/09 05:35:24 woo Exp $";
  3. #endif
  4.  
  5.  
  6. /* GNUPLOT - standard.c */
  7. /*
  8.  * Copyright (C) 1986 - 1993   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted, 
  12.  * provided that the above copyright notice appear in all copies and 
  13.  * that both that copyright notice and this permission notice appear 
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the modified code.  Modifications are to be distributed 
  18.  * as patches to released version.
  19.  *  
  20.  * This software is provided "as is" without express or implied warranty.
  21.  * 
  22.  *
  23.  * AUTHORS
  24.  * 
  25.  *   Original Software:
  26.  *     Thomas Williams,  Colin Kelley.
  27.  * 
  28.  *   Gnuplot 2.0 additions:
  29.  *       Russell Lang, Dave Kotz, John Campbell.
  30.  *
  31.  *   Gnuplot 3.0 additions:
  32.  *       Gershon Elber and many others.
  33.  * 
  34.  */
  35.  
  36. #include <math.h>
  37. #include <stdio.h>
  38. #include "plot.h"
  39.  
  40. #ifdef vms
  41. #include <errno.h>
  42. #else
  43. extern int errno;
  44. #endif /* vms */
  45.  
  46.  
  47. extern struct value stack[STACK_DEPTH];
  48. extern int s_p;
  49. extern double zero;
  50.  
  51. struct value *pop(), *Gcomplex(), *Ginteger();
  52.  
  53. double magnitude(), angle(), real(), imag();
  54.  
  55. /* The bessel function approximations here are from
  56.  * "Computer Approximations"
  57.  * by Hart, Cheney et al.
  58.  * John Wiley & Sons, 1968
  59.  */
  60.  
  61. /* There appears to be a mistake in Hart, Cheney et al. on page 149.
  62.  * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
  63.  *               Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
  64.  * In the functions below, Qn(x) is implementated using the later
  65.  * equation.
  66.  * These bessel functions are accurate to about 1e-13
  67.  */
  68.  
  69. #if defined (ATARI) && defined(__PUREC__)
  70. /* Sorry. But PUREC bugs here.
  71.  * These bessel functions are NOT accurate to about 1e-13
  72.  */
  73.  
  74. #define PI_ON_FOUR     0.785398163397448309615661
  75. #define PI_ON_TWO     1.570796326794896619231313
  76. #define THREE_PI_ON_FOUR 2.356194490192344928846982
  77. #define TWO_ON_PI     0.636619772367581343075535
  78.  
  79. static double dzero = 0.0;
  80.  
  81. /* jzero for x in [0,8]
  82.  * Index 5849, 19.22 digits precision
  83.  */
  84. static double pjzero[] = {
  85.      0.493378725179413356181681e+21,
  86.     -0.117915762910761053603844e+21,
  87.      0.638205934107235656228943e+19,
  88.     -0.136762035308817138686542e+18,
  89.      0.143435493914034611166432e+16,
  90.     -0.808522203485379387119947e+13,
  91.      0.250715828553688194555516e+11,
  92.     -0.405041237183313270636066e+8,
  93.      0.268578685698001498141585e+5
  94. };
  95.  
  96. static double qjzero[] = {
  97.      0.493378725179413356211328e+21,
  98.      0.542891838409228516020019e+19,
  99.      0.302463561670946269862733e+17,
  100.      0.112775673967979850705603e+15,
  101.      0.312304311494121317257247e+12,
  102.      0.669998767298223967181403e+9,
  103.      0.111463609846298537818240e+7,
  104.      0.136306365232897060444281e+4,
  105.      0.1e+1
  106. };
  107.  
  108. /* pzero for x in [8,inf]
  109.  * Index 6548, 18.16 digits precision
  110.  */
  111. static double ppzero[] = {
  112.      0.227790901973046843022700e+5,
  113.      0.413453866395807657967802e+5,
  114.      0.211705233808649443219340e+5,
  115.      0.348064864432492703474453e+4,
  116.      0.153762019090083542957717e+3,
  117.      0.889615484242104552360748e+0
  118. };
  119.  
  120. static double qpzero[] = {
  121.      0.227790901973046843176842e+5,
  122.      0.413704124955104166398920e+5,
  123.      0.212153505618801157304226e+5,
  124.      0.350287351382356082073561e+4,
  125.      0.157111598580808936490685e+3,
  126.      0.1e+1
  127. };
  128.  
  129. /* qzero for x in [8,inf]
  130.  * Index 6948, 18.33 digits precision
  131.  */
  132. static double pqzero[] = {
  133.     -0.892266002008000940984692e+2,
  134.     -0.185919536443429938002522e+3,
  135.     -0.111834299204827376112621e+3,
  136.     -0.223002616662141984716992e+2,
  137.     -0.124410267458356384591379e+1,
  138.     -0.8803330304868075181663e-2,
  139. };
  140.  
  141. static double qqzero[] = {
  142.      0.571050241285120619052476e+4,
  143.      0.119511315434346136469526e+5,
  144.      0.726427801692110188369134e+4,
  145.      0.148872312322837565816135e+4,
  146.      0.905937695949931258588188e+2,
  147.      0.1e+1
  148. };
  149.  
  150.  
  151. /* yzero for x in [0,8]
  152.  * Index 6245, 18.78 digits precision
  153.  */
  154. static double pyzero[] = {
  155.     -0.275028667862910958370193e+20,
  156.      0.658747327571955492599940e+20,
  157.     -0.524706558111276494129735e+19,
  158.      0.137562431639934407857134e+18,
  159.     -0.164860581718572947312208e+16,
  160.      0.102552085968639428450917e+14,
  161.     -0.343637122297904037817103e+11,
  162.      0.591521346568688965427383e+8,
  163.     -0.413703549793314855412524e+5
  164. };
  165.  
  166. static double qyzero[] = {
  167.      0.372645883898616588198998e+21,
  168.      0.419241704341083997390477e+19,
  169.      0.239288304349978185743936e+17,
  170.      0.916203803407518526248915e+14,
  171.      0.261306575504108124956848e+12,
  172.      0.579512264070072953738009e+9,
  173.      0.100170264128890626566665e+7,
  174.      0.128245277247899380417633e+4,
  175.      0.1e+1
  176. };
  177.  
  178.  
  179. /* jone for x in [0,8]
  180.  * Index 6050, 20.98 digits precision
  181.  */
  182. static double pjone[] = {
  183.      0.581199354001606143928051e+21,
  184.     -0.667210656892491629802094e+20,
  185.      0.231643358063400229793182e+19,
  186.     -0.358881756991010605074364e+17,
  187.      0.290879526383477540973760e+15,
  188.     -0.132298348033212645312547e+13,
  189.      0.341323418230170053909129e+10,
  190.     -0.469575353064299585976716e+7,
  191.      0.270112271089232341485679e+4
  192. };
  193.  
  194. static double qjone[] = {
  195.      0.116239870800321228785853e+22,
  196.      0.118577071219032099983711e+20,
  197.      0.609206139891752174610520e+17,
  198.      0.208166122130760735124018e+15,
  199.      0.524371026216764971540673e+12,
  200.      0.101386351435867398996705e+10,
  201.      0.150179359499858550592110e+7,
  202.      0.160693157348148780197092e+4,
  203.      0.1e+1
  204. };
  205.  
  206.  
  207. /* pone for x in [8,inf]
  208.  * Index 6749, 18.11 digits precision
  209.  */
  210. static double ppone[] = {
  211.      0.352246649133679798341724e+5,
  212.      0.627588452471612812690057e+5,
  213.      0.313539631109159574238670e+5,
  214.      0.498548320605943384345005e+4,
  215.      0.211152918285396238210572e+3,
  216.      0.12571716929145341558495e+1
  217. };
  218.  
  219. static double qpone[] = {
  220.      0.352246649133679798068390e+5,
  221.      0.626943469593560511888834e+5,
  222.      0.312404063819041039923016e+5,
  223.      0.493039649018108897938610e+4,
  224.      0.203077518913475932229357e+3,
  225.      0.1e+1
  226. };
  227.  
  228. /* qone for x in [8,inf]
  229.  * Index 7149, 18.28 digits precision
  230.  */
  231. static double pqone[] = {
  232.      0.351175191430355282253332e+3,
  233.      0.721039180490447503928086e+3,
  234.      0.425987301165444238988699e+3,
  235.      0.831898957673850827325226e+2,
  236.      0.45681716295512267064405e+1,
  237.      0.3532840052740123642735e-1
  238. };
  239.  
  240. static double qqone[] = {
  241.      0.749173741718091277145195e+4,
  242.      0.154141773392650970499848e+5,
  243.      0.915223170151699227059047e+4,
  244.      0.181118670055235135067242e+4,
  245.      0.103818758546213372877664e+3,
  246.      0.1e+1
  247. };
  248.  
  249.  
  250. /* yone for x in [0,8]
  251.  * Index 6444, 18.24 digits precision
  252.  */
  253. static double pyone[] = {
  254.     -0.292382196153296254310105e+20,
  255.      0.774852068218683964508809e+19,
  256.     -0.344104806308411444618546e+18,
  257.      0.591516076049007061849632e+16,
  258.     -0.486331694256717507482813e+14,
  259.      0.204969667374566218261980e+12,
  260.     -0.428947196885524880182182e+9,
  261.      0.355692400983052605669132e+6
  262. };
  263.  
  264. static double qyone[] = {
  265.      0.149131151130292035017408e+21,
  266.      0.181866284170613498688507e+19,
  267.      0.113163938269888452690508e+17,
  268.      0.475517358888813771309277e+14,
  269.      0.150022169915670898716637e+12,
  270.      0.371666079862193028559693e+9,
  271.      0.726914730719888456980191e+6,
  272.      0.107269614377892552332213e+4,
  273.      0.1e+1
  274. };
  275.  
  276. #else
  277.  
  278. #define PI_ON_FOUR       0.78539816339744830961566084581987572
  279. #define PI_ON_TWO        1.57079632679489661923131269163975144
  280. #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
  281. #define TWO_ON_PI        0.63661977236758134307553505349005744
  282.  
  283. static double dzero = 0.0;
  284.  
  285. /* jzero for x in [0,8]
  286.  * Index 5849, 19.22 digits precision
  287.  */
  288. static double pjzero[] = {
  289.      0.4933787251794133561816813446e+21,
  290.     -0.11791576291076105360384408e+21,
  291.      0.6382059341072356562289432465e+19,
  292.     -0.1367620353088171386865416609e+18,
  293.      0.1434354939140346111664316553e+16,
  294.     -0.8085222034853793871199468171e+13,
  295.      0.2507158285536881945555156435e+11,
  296.     -0.4050412371833132706360663322e+8,
  297.      0.2685786856980014981415848441e+5
  298. };
  299.  
  300. static double qjzero[] = {
  301.     0.4933787251794133562113278438e+21,
  302.     0.5428918384092285160200195092e+19,
  303.     0.3024635616709462698627330784e+17,
  304.     0.1127756739679798507056031594e+15,
  305.     0.3123043114941213172572469442e+12,
  306.     0.669998767298223967181402866e+9,
  307.     0.1114636098462985378182402543e+7,
  308.     0.1363063652328970604442810507e+4,
  309.     0.1e+1
  310. };
  311.  
  312. /* pzero for x in [8,inf]
  313.  * Index 6548, 18.16 digits precision
  314.  */
  315. static double ppzero[] = {
  316.     0.2277909019730468430227002627e+5,
  317.     0.4134538663958076579678016384e+5,
  318.     0.2117052338086494432193395727e+5,
  319.     0.348064864432492703474453111e+4,
  320.     0.15376201909008354295771715e+3,
  321.     0.889615484242104552360748e+0
  322. };
  323.  
  324. static double qpzero[] = {
  325.     0.2277909019730468431768423768e+5,
  326.     0.4137041249551041663989198384e+5,
  327.     0.2121535056188011573042256764e+5,
  328.     0.350287351382356082073561423e+4,
  329.     0.15711159858080893649068482e+3,
  330.     0.1e+1
  331. };
  332.  
  333. /* qzero for x in [8,inf]
  334.  * Index 6948, 18.33 digits precision
  335.  */
  336. static double pqzero[] = {
  337.     -0.8922660020080009409846916e+2,
  338.     -0.18591953644342993800252169e+3,
  339.     -0.11183429920482737611262123e+3,
  340.     -0.2230026166621419847169915e+2,
  341.     -0.124410267458356384591379e+1,
  342.     -0.8803330304868075181663e-2,
  343. };
  344.  
  345. static double qqzero[] = {
  346.     0.571050241285120619052476459e+4,
  347.     0.1195113154343461364695265329e+5,
  348.     0.726427801692110188369134506e+4,
  349.     0.148872312322837565816134698e+4,
  350.     0.9059376959499312585881878e+2,
  351.     0.1e+1
  352. };
  353.  
  354.  
  355. /* yzero for x in [0,8]
  356.  * Index 6245, 18.78 digits precision
  357.  */
  358. static double pyzero[] = {
  359.     -0.2750286678629109583701933175e+20,
  360.      0.6587473275719554925999402049e+20,
  361.     -0.5247065581112764941297350814e+19,
  362.      0.1375624316399344078571335453e+18,
  363.     -0.1648605817185729473122082537e+16,
  364.      0.1025520859686394284509167421e+14,
  365.     -0.3436371222979040378171030138e+11,
  366.      0.5915213465686889654273830069e+8,
  367.     -0.4137035497933148554125235152e+5
  368. };
  369.  
  370. static double qyzero[] = {
  371.     0.3726458838986165881989980739e+21,
  372.     0.4192417043410839973904769661e+19,
  373.     0.2392883043499781857439356652e+17,
  374.     0.9162038034075185262489147968e+14,
  375.     0.2613065755041081249568482092e+12,
  376.     0.5795122640700729537380087915e+9,
  377.     0.1001702641288906265666651753e+7,
  378.     0.1282452772478993804176329391e+4,
  379.     0.1e+1
  380. };
  381.  
  382.  
  383. /* jone for x in [0,8]
  384.  * Index 6050, 20.98 digits precision
  385.  */
  386. static double pjone[] = {
  387.      0.581199354001606143928050809e+21,
  388.     -0.6672106568924916298020941484e+20,
  389.      0.2316433580634002297931815435e+19,
  390.     -0.3588817569910106050743641413e+17,
  391.      0.2908795263834775409737601689e+15,
  392.     -0.1322983480332126453125473247e+13,
  393.      0.3413234182301700539091292655e+10,
  394.     -0.4695753530642995859767162166e+7,
  395.      0.270112271089232341485679099e+4
  396. };
  397.  
  398. static double qjone[] = {
  399.     0.11623987080032122878585294e+22,
  400.     0.1185770712190320999837113348e+20,
  401.     0.6092061398917521746105196863e+17,
  402.     0.2081661221307607351240184229e+15,
  403.     0.5243710262167649715406728642e+12,
  404.     0.1013863514358673989967045588e+10,
  405.     0.1501793594998585505921097578e+7,
  406.     0.1606931573481487801970916749e+4,
  407.     0.1e+1
  408. };
  409.  
  410.  
  411. /* pone for x in [8,inf]
  412.  * Index 6749, 18.11 digits precision
  413.  */
  414. static double ppone[] = {
  415.     0.352246649133679798341724373e+5,
  416.     0.62758845247161281269005675e+5,
  417.     0.313539631109159574238669888e+5,
  418.     0.49854832060594338434500455e+4,
  419.     0.2111529182853962382105718e+3,
  420.     0.12571716929145341558495e+1
  421. };
  422.  
  423. static double qpone[] = {
  424.     0.352246649133679798068390431e+5,
  425.     0.626943469593560511888833731e+5,
  426.     0.312404063819041039923015703e+5,
  427.     0.4930396490181088979386097e+4,
  428.     0.2030775189134759322293574e+3,
  429.     0.1e+1
  430. };
  431.  
  432. /* qone for x in [8,inf]
  433.  * Index 7149, 18.28 digits precision
  434.  */
  435. static double pqone[] = {
  436.     0.3511751914303552822533318e+3,
  437.     0.7210391804904475039280863e+3,
  438.     0.4259873011654442389886993e+3,
  439.     0.831898957673850827325226e+2,
  440.     0.45681716295512267064405e+1,
  441.     0.3532840052740123642735e-1
  442. };
  443.  
  444. static double qqone[] = {
  445.     0.74917374171809127714519505e+4,
  446.     0.154141773392650970499848051e+5,
  447.     0.91522317015169922705904727e+4,
  448.     0.18111867005523513506724158e+4,
  449.     0.1038187585462133728776636e+3,
  450.     0.1e+1
  451. };
  452.  
  453.  
  454. /* yone for x in [0,8]
  455.  * Index 6444, 18.24 digits precision
  456.  */
  457. static double pyone[] = {
  458.     -0.2923821961532962543101048748e+20,
  459.      0.7748520682186839645088094202e+19,
  460.     -0.3441048063084114446185461344e+18,
  461.      0.5915160760490070618496315281e+16,
  462.     -0.4863316942567175074828129117e+14,
  463.      0.2049696673745662182619800495e+12,
  464.     -0.4289471968855248801821819588e+9,
  465.      0.3556924009830526056691325215e+6
  466. };
  467.  
  468. static double qyone[] = {
  469.     0.1491311511302920350174081355e+21,
  470.     0.1818662841706134986885065935e+19,
  471.     0.113163938269888452690508283e+17,
  472.     0.4755173588888137713092774006e+14,
  473.     0.1500221699156708987166369115e+12,
  474.     0.3716660798621930285596927703e+9,
  475.     0.726914730719888456980191315e+6,
  476.     0.10726961437789255233221267e+4,
  477.     0.1e+1
  478. };
  479.  
  480. #endif /* ATARI && __PUREC__ */
  481.  
  482. f_real()
  483. {
  484. struct value a;
  485.     push( Gcomplex(&a,real(pop(&a)), 0.0) );
  486. }
  487.  
  488. f_imag()
  489. {
  490. struct value a;
  491.     push( Gcomplex(&a,imag(pop(&a)), 0.0) );
  492. }
  493.  
  494. f_arg()
  495. {
  496. struct value a;
  497.     push( Gcomplex(&a,angle(pop(&a)), 0.0) );
  498. }
  499.  
  500. f_conjg()
  501. {
  502. struct value a;
  503.     (void) pop(&a);
  504.     push( Gcomplex(&a,real(&a),-imag(&a) ));
  505. }
  506.  
  507. f_sin()
  508. {
  509. struct value a;
  510.     (void) pop(&a);
  511.     push( Gcomplex(&a,sin(real(&a))*cosh(imag(&a)), cos(real(&a))*sinh(imag(&a))) );
  512. }
  513.  
  514. f_cos()
  515. {
  516. struct value a;
  517.     (void) pop(&a);
  518.     push( Gcomplex(&a,cos(real(&a))*cosh(imag(&a)), -sin(real(&a))*sinh(imag(&a))));
  519. }
  520.  
  521. f_tan()
  522. {
  523. struct value a;
  524. register double den;
  525.     (void) pop(&a);
  526.     if (imag(&a) == 0.0)
  527.         push( Gcomplex(&a,tan(real(&a)),0.0) );
  528.     else {
  529.         den = cos(2*real(&a))+cosh(2*imag(&a));
  530.         if (den == 0.0) {
  531.             undefined = TRUE;
  532.             push( &a );
  533.         }
  534.         else
  535.             push( Gcomplex(&a,sin(2*real(&a))/den, sinh(2*imag(&a))/den) );
  536.     }
  537. }
  538.  
  539. f_asin()
  540. {
  541. struct value a;
  542. register double alpha, beta, x, y;
  543.     (void) pop(&a);
  544.     x = real(&a); y = imag(&a);
  545.     if (y == 0.0) {
  546.         if (fabs(x) > 1.0) {
  547.             undefined = TRUE;
  548.             push(Gcomplex(&a,0.0, 0.0));
  549.         } else
  550.             push( Gcomplex(&a,asin(x),0.0) );
  551.     } else {
  552.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  553.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  554.         push( Gcomplex(&a,asin(beta), log(alpha + sqrt(alpha*alpha-1))) );
  555.     }
  556. }
  557.  
  558. f_acos()
  559. {
  560. struct value a;
  561. register double alpha, beta, x, y;
  562.     (void) pop(&a);
  563.     x = real(&a); y = imag(&a);
  564.     if (y == 0.0) {
  565.         if (fabs(x) > 1.0) {
  566.             undefined = TRUE;
  567.             push(Gcomplex(&a,0.0, 0.0));
  568.         } else
  569.             push( Gcomplex(&a,acos(x),0.0) );
  570.     } else {
  571.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  572.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  573.         push( Gcomplex(&a,acos(beta), log(alpha + sqrt(alpha*alpha-1))) );
  574.     }
  575. }
  576.  
  577. f_atan()
  578. {
  579. struct value a;
  580. register double x, y, u, v, w, z;
  581.     (void) pop(&a);
  582.     x = real(&a); y = imag(&a);
  583.     if (y == 0.0)
  584.         push( Gcomplex(&a,atan(x), 0.0) );
  585.     else if (x == 0.0 && fabs(y) == 1.0) {
  586.         undefined = TRUE;
  587.         push(Gcomplex(&a,0.0, 0.0));
  588.     } else {
  589.             if (x >= 0) {
  590.                 u = x;
  591.             v = y;
  592.         } else {
  593.                 u = -x;
  594.             v = -y;
  595.         }
  596.         
  597.             z = atan(2*u/(1-u*u-v*v));
  598.         w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
  599.         if (z < 0)
  600.                 z = z + 2*PI_ON_TWO;
  601.         if (x < 0) {
  602.                 z = -z;
  603.             w = -w;
  604.         }
  605.         push( Gcomplex(&a,0.5*z, w) );
  606.     }
  607. }
  608.  
  609. f_sinh()
  610. {
  611. struct value a;
  612.     (void) pop(&a);
  613.     push( Gcomplex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
  614. }
  615.  
  616. f_cosh()
  617. {
  618. struct value a;
  619.     (void) pop(&a);
  620.     push( Gcomplex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
  621. }
  622.  
  623. f_tanh()
  624. {
  625. struct value a;
  626. register double den;
  627.     (void) pop(&a);
  628.     den = cosh(2*real(&a)) + cos(2*imag(&a));
  629.     push( Gcomplex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
  630. }
  631.  
  632. f_int()
  633. {
  634. struct value a;
  635.     push( Ginteger(&a,(int)real(pop(&a))) );
  636. }
  637.  
  638.  
  639. f_abs()
  640. {
  641. struct value a;
  642.     (void) pop(&a);
  643.     switch (a.type) {
  644.         case INTGR:
  645.             push( Ginteger(&a,abs(a.v.int_val)) );            
  646.             break;
  647.         case CMPLX:
  648.             push( Gcomplex(&a,magnitude(&a), 0.0) );
  649.     }
  650. }
  651.  
  652. f_sgn()
  653. {
  654. struct value a;
  655.     (void) pop(&a);
  656.     switch(a.type) {
  657.         case INTGR:
  658.             push( Ginteger(&a,(a.v.int_val > 0) ? 1 : 
  659.                     (a.v.int_val < 0) ? -1 : 0) );
  660.             break;
  661.         case CMPLX:
  662.             push( Ginteger(&a,(a.v.cmplx_val.real > 0.0) ? 1 : 
  663.                     (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
  664.             break;
  665.     }
  666. }
  667.  
  668.  
  669. f_sqrt()
  670. {
  671. struct value a;
  672. register double mag, ang;
  673.     (void) pop(&a);
  674.     mag = sqrt(magnitude(&a));
  675.     if (imag(&a) == 0.0 && real(&a) < 0.0)
  676.         push( Gcomplex(&a,0.0,mag) );
  677.     else
  678.     {
  679.         if ( (ang = angle(&a)) < 0.0)
  680.             ang += 2*Pi;
  681.         ang /= 2;
  682.         push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
  683.     }
  684. }
  685.  
  686.  
  687. f_exp()
  688. {
  689. struct value a;
  690. register double mag, ang;
  691.     (void) pop(&a);
  692.     mag = exp(real(&a));
  693.     ang = imag(&a);
  694.     push( Gcomplex(&a,mag*cos(ang), mag*sin(ang)) );
  695. }
  696.  
  697.  
  698. f_log10()
  699. {
  700. struct value a;
  701. register double l10;;
  702.     (void) pop(&a);
  703.     l10 = log(10.0);    /***** replace with a constant! ******/
  704.     push( Gcomplex(&a,log(magnitude(&a))/l10, angle(&a)/l10) );
  705. }
  706.  
  707.  
  708. f_log()
  709. {
  710. struct value a;
  711.     (void) pop(&a);
  712.     push( Gcomplex(&a,log(magnitude(&a)), angle(&a)) );
  713. }
  714.  
  715.  
  716. f_floor()
  717. {
  718. struct value a;
  719.  
  720.     (void) pop(&a);
  721.     switch (a.type) {
  722.         case INTGR:
  723.             push( Ginteger(&a,(int)floor((double)a.v.int_val)));            
  724.             break;
  725.         case CMPLX:
  726.             push( Ginteger(&a,(int)floor(a.v.cmplx_val.real)));
  727.     }
  728. }
  729.  
  730.  
  731. f_ceil()
  732. {
  733. struct value a;
  734.  
  735.     (void) pop(&a);
  736.     switch (a.type) {
  737.         case INTGR:
  738.             push( Ginteger(&a,(int)ceil((double)a.v.int_val)));            
  739.             break;
  740.         case CMPLX:
  741.             push( Ginteger(&a,(int)ceil(a.v.cmplx_val.real)));
  742.     }
  743. }
  744.  
  745. /* bessel function approximations */
  746. double jzero(x)
  747. double x;
  748. {
  749. double p, q, x2;
  750. int n;
  751.  
  752.     x2 = x * x;
  753.     p = pjzero[8];
  754.     q = qjzero[8];
  755.     for (n=7; n>=0; n--) {
  756.         p = p*x2 + pjzero[n];
  757.         q = q*x2 + qjzero[n];
  758.     }
  759.     return(p/q);
  760. }
  761.  
  762. double pzero(x)
  763. double x;
  764. {
  765. double p, q, z, z2;
  766. int n;
  767.  
  768.     z = 8.0 / x;
  769.     z2 = z * z;
  770.     p = ppzero[5];
  771.     q = qpzero[5];
  772.     for (n=4; n>=0; n--) {
  773.         p = p*z2 + ppzero[n];
  774.         q = q*z2 + qpzero[n];
  775.     }
  776.     return(p/q);
  777. }
  778.  
  779. double qzero(x)
  780. double x;
  781. {
  782. double p, q, z, z2;
  783. int n;
  784.  
  785.     z = 8.0 / x;
  786.     z2 = z * z;
  787.     p = pqzero[5];
  788.     q = qqzero[5];
  789.     for (n=4; n>=0; n--) {
  790.         p = p*z2 + pqzero[n];
  791.         q = q*z2 + qqzero[n];
  792.     }
  793.     return(p/q);
  794. }
  795.  
  796. double yzero(x)
  797. double x;
  798. {
  799. double p, q, x2;
  800. int n;
  801.  
  802.     x2 = x * x;
  803.     p = pyzero[8];
  804.     q = qyzero[8];
  805.     for (n=7; n>=0; n--) {
  806.         p = p*x2 + pyzero[n];
  807.         q = q*x2 + qyzero[n];
  808.     }
  809.     return(p/q);
  810. }
  811.  
  812. double rj0(x)
  813. double x;
  814. {
  815.     if ( x <= 0.0 )
  816.         x = -x;
  817.     if ( x < 8.0 )
  818.         return(jzero(x));
  819.     else
  820.         return( sqrt(TWO_ON_PI/x) *
  821.             (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
  822.  
  823. }
  824.  
  825. double ry0(x)
  826. double x;
  827. {
  828.     if ( x < 0.0 )
  829.         return(dzero/dzero); /* error */
  830.     if ( x < 8.0 )
  831.         return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
  832.     else
  833.         return( sqrt(TWO_ON_PI/x) *
  834.             (pzero(x)*sin(x-PI_ON_FOUR) + 
  835.             (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
  836.  
  837. }
  838.  
  839.  
  840. double jone(x)
  841. double x;
  842. {
  843. double p, q, x2;
  844. int n;
  845.  
  846.     x2 = x * x;
  847.     p = pjone[8];
  848.     q = qjone[8];
  849.     for (n=7; n>=0; n--) {
  850.         p = p*x2 + pjone[n];
  851.         q = q*x2 + qjone[n];
  852.     }
  853.     return(p/q);
  854. }
  855.  
  856. double pone(x)
  857. double x;
  858. {
  859. double p, q, z, z2;
  860. int n;
  861.  
  862.     z = 8.0 / x;
  863.     z2 = z * z;
  864.     p = ppone[5];
  865.     q = qpone[5];
  866.     for (n=4; n>=0; n--) {
  867.         p = p*z2 + ppone[n];
  868.         q = q*z2 + qpone[n];
  869.     }
  870.     return(p/q);
  871. }
  872.  
  873. double qone(x)
  874. double x;
  875. {
  876. double p, q, z, z2;
  877. int n;
  878.  
  879.     z = 8.0 / x;
  880.     z2 = z * z;
  881.     p = pqone[5];
  882.     q = qqone[5];
  883.     for (n=4; n>=0; n--) {
  884.         p = p*z2 + pqone[n];
  885.         q = q*z2 + qqone[n];
  886.     }
  887.     return(p/q);
  888. }
  889.  
  890. double yone(x)
  891. double x;
  892. {
  893. double p, q, x2;
  894. int n;
  895.  
  896.     x2 = x * x;
  897.     p = 0.0;
  898.     q = qyone[8];
  899.     for (n=7; n>=0; n--) {
  900.         p = p*x2 + pyone[n];
  901.         q = q*x2 + qyone[n];
  902.     }
  903.     return(p/q);
  904. }
  905.  
  906. double rj1(x)
  907. double x;
  908. {
  909. double v,w;
  910.     v = x;
  911.     if ( x < 0.0 )
  912.         x = -x;
  913.     if ( x < 8.0 )
  914.         return(v*jone(x));
  915.     else {
  916.         w = sqrt(TWO_ON_PI/x) *
  917.             (pone(x)*cos(x-THREE_PI_ON_FOUR) - 
  918.                8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
  919.         if (v < 0.0)
  920.             w = -w;
  921.         return( w );
  922.     }
  923. }
  924.  
  925. double ry1(x)
  926. double x;
  927. {
  928.     if ( x <= 0.0 )
  929.         return(dzero/dzero); /* error */
  930.     if ( x < 8.0 )
  931.         return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
  932.     else
  933.         return( sqrt(TWO_ON_PI/x) *
  934.             (pone(x)*sin(x-THREE_PI_ON_FOUR) + 
  935.             (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
  936. }
  937.  
  938.  
  939. f_besj0()    
  940. {
  941. struct value a;
  942.     (void) pop(&a);
  943.     if (fabs(imag(&a)) > zero)
  944.         int_error("can only do bessel functions of reals",NO_CARET);
  945.     push( Gcomplex(&a,rj0(real(&a)),0.0) );
  946. }
  947.  
  948.  
  949. f_besj1()    
  950. {
  951. struct value a;
  952.     (void) pop(&a);
  953.     if (fabs(imag(&a)) > zero)
  954.         int_error("can only do bessel functions of reals",NO_CARET);
  955.     push( Gcomplex(&a,rj1(real(&a)),0.0) );
  956. }
  957.  
  958.  
  959. f_besy0()    
  960. {
  961. struct value a;
  962.     (void) pop(&a);
  963.     if (fabs(imag(&a)) > zero)
  964.         int_error("can only do bessel functions of reals",NO_CARET);
  965.     if (real(&a) > 0.0)
  966.         push( Gcomplex(&a,ry0(real(&a)),0.0) );
  967.     else {
  968.         push( Gcomplex(&a,0.0,0.0) );
  969.         undefined = TRUE ;
  970.     }
  971. }
  972.  
  973.  
  974. f_besy1()    
  975. {
  976. struct value a;
  977.     (void) pop(&a);
  978.     if (fabs(imag(&a)) > zero)
  979.         int_error("can only do bessel functions of reals",NO_CARET);
  980.     if (real(&a) > 0.0)
  981.         push( Gcomplex(&a,ry1(real(&a)),0.0) );
  982.     else {
  983.         push( Gcomplex(&a,0.0,0.0) );
  984.         undefined = TRUE ;
  985.     }
  986. }
  987.